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We propose a new type of the Ensemble Kalman Filter (EnKF), which uses 
the Fast Fourier Transform (FFT) for covariance estimation from a very small 
^ ensemble with automatic tapering, and for a fast computation of the analy- 

sis ensemble by convolution, avoiding the need to solve a. sparse SA'stem with 
the tapered matrix. The method is combined with the niorphiug EnKF to 
enable the correction of position errors, in addition to amplitude errors, and 
demonstrated on WRF-Fire, the Weather Research Forecasting (WRF) model 
coupled with a fire spread model implemented by the level set method. 
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^ 1. Introduction 

00 

00 Data assimilation is a statistical technique to modify the state of a running 

model in response to data, based on sequential Bayesian estimation.^ The 
EnKF^ accesses the model only as a black box. It is suitable for a wide 
range of problems and relatively easy to combine with existing modeling 
software. However, the ensemble size needed can be large, and the amount 

^-H of computation required in EnKF can be significant when the EnKF is 

^ modified to suppress spurious long-range correlations. We propose a new 

: variant of EnKF that can overcome these difficiilties by the use of FFT. The 

rS new method is applied to the morphing EnKF,"^ which is needed for the fire 

application, but the FFT EnKF can be used with or without morphing. 



2. Data Assimilation 

We first briefly describe the EnKF for reference. The EnKF advances in time 

an ensemble of simulations ui , . . . , uj^ , which approximates the probability 
distribution of the model state u. The simulations are then advanced in 
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time until the analysis time, when new data d arrives. It is assumed that the 
data error is known, d — Hu ~ N (0, R) given u, where H is the observation 
operator and R is the data error covariance matrix. The ensemble, now 
called the forecast ensemble, is combined with the data to give the analysis 
ensemble by the EnKF formulas^ 



where Cat is an estimate of the covariance of the model state, and Ck is a 
random perturbation ~ N{0,R). The analysis ensemble 
then used as the initial condition for a new simulation, advanced to a new 
analysis time, and the analysis cycle repeats. 

In the standard EnKF,^ Cat is the sample covariance computed from the 
ensemble. Under suitable assumptions, it can be proved that the ensemble 
converges for large A?' to a sample from the Kalman filtering distribution 
and Cjv converges to the true covariance.^ 

2.1. FFT EnKF and Covariance Estimation by FFT 

The EnKF formulation (1) relies on linear algebra only and it is oblivious 
to any structure of the model state. Large ensembles (30-100 and more) 
are often needed.^ Wc are interested in obtaining a reasonable approxima- 
tion of the covariance matrix from a very small sample. For this, we take 
advantage of the fact that the simulation state w is a block vector, where 
the blocks are values of the modelled physical fields on grids of points in 
a spatial domain fi, which are (discrete versions of) smooth random func- 
tions, i.e., realizations of random fields. The covariance of a random field 
drops off quickly with distance, and it is often the case in geostatistics that 
random fields are stationary, that is, the covariance between two points de- 
pends on their distance vector only.^ But for a small sample, the ensemble 
covariance is a matrix of low rank with large off-diagonal elements even 
at a great distance. Therefore, localization techniques are used, such as 
tapering, which consists of multiplication of the terms of the sample covari- 
ance by a fixed function to force the dropoff of the covariance away from 
the diagonal, resulting in a more accurate approximation of covariance for 
small samples.'' However, solving a system with the resulting approximate 
covariance matrix is expensive, because efficient dense linear algebra, which 
relies on the representation of the sample covariance matrix as the product 
of two rectangular matrices,^ can no longer be used. 

For simplicity, we explain the FFT EnKF in the ID case. Higher- 
dimensional cases work exactly the same. Also, we consider first the case 




(1) 
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when the model state consists of one block only. 

The basic operation in the EnKF (1) is the multiplication v = Cnu of 
a vector u by an approximate covariance matrix Cjv- Denote by u{xi) the 
entry of vector w, corresponding to node Xi and suppose for the moment 
that the random field is stationary, that is, its covariance matrix satisfies 
C {xi,Xj) = c{xi — Xj) for some covariance function c. Then v is the con- 
volution 

V {Xi) = C (Xi, Xj) U (Xj) = U (Xj) C {Xi — Xj) . 

j 3 

The discrete Fourier transform changes convolution to multiplication en- 
try by entry. Hence, if the random field is stationary, multiplication by 
its covariance matrix becomes multiplication by a diagonal matrix in the 
frequency domain. 

The proposed method of covariance estimation consists of computing the 
sample covariance of the ensemble in the frequency domain, and neglecting 
all of its off-diagonal terms. This is justified by the assumption that the 
covariance depends mainly on the distance. 

(1) Given ensemble [itfe], apply a FFT operator F to each member to obtain 
the member Uk = Fu^ in the frequency domain. 

(2) Compute the approximate forecast covariance matrix Cm the fre- 
quency domain as the diagonal matrix with the diagonal entries ci 
equal to the diagonal entries of the sample covariance of the ensemble 



Ci 



N 



fc=l k=l 



^E"^^- (2) 



Multiplication by the approximate covariance matrix Cjv then becomes 
in the frequency domain 

u = Cnv ■^=^ u = Fu, v = c»u, V = F~^v, 

where • is cntry-by-entry multiplication, (c • u)^ = CiUi. In the important 
case H = I and R = rl (the whole state is observed and the data errors 
are uncorrelated and the same at every point), considered here, the EnKF 
(1) in the frequency domain reduces to entry- by-entry operations, 

ul=Uk + c • {c + r)~^ • {d + Cfe - m{ j . (3) 
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In an application, the state has multiple variables. The state vector, its 
covariance, and the observation matrix then have the block form 





■yd)- 






. ^(IM)- 




u = 




, c = 


C(M1) . 




, H = [if (1) • • • ifW] . (4) 



Assume that the first variable is observed, then H^^^ = I, H^"^^ = 0,. . . , 
jj(M) ^ 0. The EnKF (1) then simplifies to 

4^')'" = 4^'^ + C^'^ {C^^^^ + r) (d + eu- «W) , i = 1, . . . , M, (5) 

which becomes in the frequency domain 

4)-'^ = 4) + c<^i) . (c<") +ry\(d + eu- Uk) , (6) 

where the spectral cross-covariance between field j and field 1 is approxi- 
mated from 

fc=i fc=i 

2.2. Morphing EnKF 

To treat position errors in addition to amplitude errors, FFT EnKF is com- 
bined with the morphing EnKF.'^'* The method uses an additional ensemble 
member ujv+i, called the reference member. Given an initial state u, the 
initial ensemble is given by Ujv+i = u and 

4^ = (4li + r-r) o{i+n),k = i,...,N, (8) 

where r^^ are random smooth functions on fi, are random smooth map- 
pings Th : Cl ^ SI, and o denotes composition. Thus, the initial ensemble 

has both amplitude and position variability, and the position change is the 
same for all blocks. Random smooth functions and mappings are generated 
by FFT as a Fourier series with random coefficients that decay quickly with 
frequency. 

The data d is an observation of w'^^^ and it is expected that it differs 
from the model in amplitude as well as in the position of significant features, 
such as firelines. The first blocks of ui , . . . , ujv and d are then registered 
against the first block of the reference member utv+i. We find registration 
mappings : — > $7, fc = 0, . . . , such that 

4'^ ~"^li°U + ^fe)' Tk^O, VTk^O, k = 0,...,N, 
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where d = Define the registration residuals r^^ = u^i^^ o (7 + T/j)~^ — 
u^N+i^ k = 0, . . . , N. The m,orphing transform maps each ensemble member 
Wfeinto the extended state vector 

Uk^Uk = M„„+, (ufe) = (Tfe, ri'\ . . . , ri^^) . (9) 

Similarly, the data becomes the extended data vector d = {Tq^t'i^^^. 
The FFT EnKF method (6) is applied to the transformed ensemble 
«i, . . . , UN with the observation operator given by(T, r'^^\ . . . , r^^^) i-)- 
(T, r^^'). The cross-covariances between x and y components of T and r'^^ 
are neglected, so the covariance C^^^-* in (5) consists of three diagonal matri- 
ces, and (6) applies. The new transformed reference member is obtained as 
uj^_^_i = jfJ2k=i''^k ^^'^ analysis ensemble ui, . . . ,un+i by the inverse 
morphing transform 

A^««V(«fe) = (4Vi+r,"'''^)°(/ + T,^), fc = l,...,iV+l. (10) 



3. The Wildland Fire Model 

We only summarize the model very briefly. See Ref. 8 for further details, 
references, and acknowledgements. 

The fire model runs on a 2D rectangular grid on the Earth surface, called 
the fire grid. The model postulates fire line propagation speed as a function 
of wind and terrain slope, and an exponential decay of fuel by combustion 
after the ignition. The fire area is represented by a level set function as the 
set of points where the level set function is negative. The level set function 
satisfies a partial differential equation, which is solved numerically on the 
fire grid by explicit time stepping. Ignition is achieved by setting the value 
of the level set function so that it is negative in the ignition region. The 
state of the fire model consists of the level set function and the ignition 
time at the fire grid nodes. Other derived quantities included in the state 
vector are the fuel fraction burned and the heat flux from the combustion. 

The fire model is coupled with the Weather Research Forecasting 
(WRF) atmospheric model, ^ which runs on a much coarser 3D grid. In 
every time step, the fire model inputs the horizontal wind velocity, and 
the heat and vapor fluxes output from the fire model are inserted into the 
atmospheric model. The fire model is distributed as WRF-Fire with WRF, 
and the latest development version is also available from the authors. 
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4. Computational Results 

Wc have used the optimization method from Ref. 3 for registration. We 
have used the real sin FFT, which forces zero change on the boundary. 
We have used a standard ideal problem distributed with WRF-Fire. The 
model has a 420 x 420 fire mesh and a 42 x 42 x 41 atmospheric mesh. The 
fuel was the same on the whole domain. The model was initialized with 
the wind blowing diagonally across the grid, and two line ignitions and one 
circle ignition occur within the first 4 seconds of simulation time. After 
one minute of simulation time, when the fire was established and one of 
the line fires has merged with the circular fire, the simulation was stopped 
and an initial ensemble was generated by random smooth perturbation 
both in position and in amplitude. Artificial data was created by a similar 
perturbation. The forecast was taken the same as the initial ensemble. The 
described data assimilation algorithm was then applied with 5 members, 
with the results shown in Fig. 1 for the morphing EnKF and Fig. 2 for the 
morphing FFT EnKF. Wc sec that the EnKF was not able to approach 
that data at all with such a small ensemble, while the FFT EnKF delivered 
an ensemble around the correct data shape. 

5. Conclusion 

We have shown that the morphing FFT EnKF is capable of data assimila- 
tion in a wildfire simulation, which exhibits sharp boundaries and coherent 
features. We have shown that the FFT EnKF can deliver acceptable re- 
sults with a very small ensemble (5 members), unlike the standard EnKF, 
which is known to work with morphing for this application, but only with 
a much larger ensemble.^ Further development, involving multiple variables 
and multiple analysis cycles, will be presented elsewhere. 
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(c) Analysis member 3 
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Fig. 1. The morphing EnKF with 5 ensemble members, applied to the ground heat flux 
from the WRF-Fire model. The ensemble size is not sufficient, the correct analysis is not 
even approximately in the span of the forecast, and the EnKF cannot reach it. 
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Fig. 2. The morphing FFT EiiKF with 5 ensemble members, appUed to the ground 
heat flux from the WRF-Fire modeh The analysis ensemble moved towards the data. 



